Quantum phase transition in a two-dimensional system of dipoles 
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The ground-state phase diagram of a two-dimensional Bose system with dipole-dipole interactions 
is studied by means of quantum Monte Carlo technique. Our calculation predicts a quantum phase 
transition from gas to solid phase when the density increases. In the gas phase the condensate 
fraction is calculated as a function of the density. Using Feynman approximation, the collective 
excitation branch is studied and appearance of a roton minimum is observed. Results of the static 
structure factor at both sides of the gas-solid phase are also presented. The Lindeman ratio at the 
transition point comes to be 7 = 0.230(6). The condensate fraction in the gas phase is estimated as 
a function of the density. 



The chromium atom has exceptionally large permanent 
dipole moment and recent realization of Bose-Einstein 
condensation of 52 Cr 1] has stimulated great interest in 
properties of dipolar systems at low temperatures. It 
was observed that dipolar forces lead to anisotropic 
deformation during expansion of the condensate. In the 
experiments [l], [2] , the dipolar forces were competing with 
short-range scattering. The latter, in principle, can be 
removed by tuning the s-wave scattering length to zero 
by Feshbach resonance [3(- This would lead to an essen- 
tially pure system of dipoles. On the other hand, low- 
dimensional systems can be realized by making the con- 
finement in one (or two) directions so tight, that no exci- 
tations of the levels of the tight confinement are possible 
and the system is dynamically two- (or one-) dimensional. 

A major development has also been done in the present 
years towards the realization of excitons at temperatures 
close to the Bose-Einstein condensation temperature^. 
An exciton is much more stable if the electron is spa- 
tially separated from the hole (spatially separated exci- 
tons). Such an exciton can be modeled as a dipole. If 
the excitons are in two coupled quantum wells they can 
be treated effectively as two-dimensional if the size of an 
exciton is small compared to the mean distance between 
excitons. 

One might expect to find a phase transition from gas 
phase to a crystal one at large density. As the conden- 
sate fraction is small at the transition point, perturbative 
theories, like Gross-Pitaevskii [f| or Bogoliubov 0,0] ap- 
proaches, will fail to describe accurately this transition. 
One has to use ab initio numerical methods to address 
this quantum many-body problem. Recently a trapped 
system of two-dimensional dipoles has been studied by 
Path Integral Monte Carlo method0 and mesoscopic 
analog of crystallization has been found. Trapped dipoles 
with s-wave scattering were investigated^]. So far, there 
have been no full quantum microscopic computations of 
the properties of a homogeneous system of dipoles. 

The Hamiltonian of a homogeneous system of N 



bosonic dipoles has the form 
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where M is the dipole mass and n, i = 1,N are the po- 
sitions of dipoles. The expression for the coupling con- 
stant Cdd depends on the nature of the dipole-dipole in- 
teraction. Two possible physical realizations of a two- 
dimensional system of dipoles can be considered: 

1) Cold atoms with permanent dipole moments m 
aligned perpendicularly to the plane of confinement by 
an external magnetic field. In this case Cdd = Mo m2 > 
where [1q is the permeability of free space. Alternatively, 
the electric dipole moment can be induced by an electric 
field E, then the coupling constant is Cdd = E 2 a 2 /eo, 
where a is the static polarizability and eo the permittiv- 
ity of free space. 

2) Spatially indirect excitons. In two coupled quantum 
wells, one containing only holes, and the other only elec- 
trons, holes and electrons might couple forming indirect 
excitons. Another possible realization is a single quan- 
tum well where the exciton dipole moment is induced by 
normal electric field. Spatial separation between hole and 
electron suppress recombination and greatly increases the 
lifetime of an exciton. If the separation between excitons 
are greater than the electron - hole separation D, indirect 
excitons can be approximated as bosons with with dipo- 
lar moment oriented perpendicularly to the plane. In this 
case Cdd = e 2 D 2 /e, where e is an electron's charge and 
e is the dielectric constant of the semiconductor. Study 
of 2D indirect exciton systems in fact is a hot topic and 
this problem has bej3n addressed both theoretically (hJ 
and experimentally 

The Hamiltonian |T]) can be written in dimcnsionless 
form by expressing all distances in units of characteris- 
tic length ro = MCdd/i^nh 2 ) and energies in units of 
£0 = ti 2 /Mr 2 .. Properties of a homogeneous system are 
governed by the dimcnsionless parameter nr 2 . (dimcn- 
sionless density) with n being density of the system. 
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In this paper we present a complete study of the 
phase diagram of two-dimensional bosonic particles with 
dipole-dipole interactions at zero temperature. We re- 
sort to the Diffusion Monte Carlo (DMC) method[ll| 
in order to find the ground-state energy and correlation 
functions of the many-body Hamiltonian (JXJ) . Within 
DMC method the Schrodinger equation is solved in 
imaginary time for the product of the ground-state 
wave function and a trial (or guiding) wave function, 
which we chose of Bijl-Jastrow form: ^t{^i, ■ tn) 
1 1 , • /i(ri) n^fe Mr; - r k |). In the gas phase the den- 
sity profile is constant and /i(r) = const. When two 
particles closely approach, the influence of other parti- 
cles can be neglected and we expect that /2(f) is well 
approximated by the solution of the two-body scattering 
problem. We choose the short distance part of the two- 
body correlation term as f%(r) = CiK (2/y/r), where 
Kq(t) is modified Bessel function of the second kind 
and C'i is some constant. At large distances, instead, 
contribution from other particles must not be neglected 
and collective behavior (phonons) is expected. From 
hydrodynamic theory it was shown [l2| that in a two- 
dimensional system the phononic long-range part of the 
wave function decays as /2(f) oc exp(— const/r). We 
use this functional dependence in a symmetrized form 
f%(r) = Ci exp(— 63/r) exp(— C^/(L~ r)), which ensures 
that = 0. We match the short- and long- range 

parts of fa(r) and demand continuity of the function 
and its first derivative at a matching distance R ma tch- 
This, together with the condition / 2 (L/2) = 1, fixes con- 
stants 61,62,63. In the solid phase we use the same 
/2(f) and add the one-body Gaussian localization term 
fi( r i) = exp[— a(ri — rf r ) 2 ],i — l,N, where a is the lo- 
calization strength and rf position of the lattice site. 
Variational parameters R m atch and a are chosen by min- 
imizing the variational energy. 

We model the homogeneous system at density n by 
considering N particles in a simulation box with periodic 
boundary conditions 
chosen in such a way that n 



E/N 



The size of the box L x x L y is 
N/(L x L y ). In the gas 



phase a square box L x — L y is considered, while in the 
solid phase we ensure that each of the box sides is a 
multiple of an elementary cell size in a triangular lattice. 

The advantage of DMC method is that it gives essen- 
tially exact energy within some statistical uncertainties. 
In a system of dipoles special attention should be paid 
to an appropriate extrapolation to the infinite system. 
The reason for that is that while 1/r 3 dipole-dipole in- 
teraction is not yet a long-range one in a two-dimensional 
system [3], the long-range decay is already quite slow. 
The finite-size effects in energy can be significantly re- 
duced by adding the tail energy 



E tail (n,L) 
N 



V (r)g2(r)2nrdr 



(2) 
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FIG. 1: (color online) Energy per particle of the dipole system 
as a function of the dimensionless density nr 2 ,: (blue) solid 
squares - gas phase, (red) open circles - solid phase. Energy 
is measured in units [h 2 / (Mr^)] / (m , o) 3//2 . Inset: an example 
of the finite-size dependence for the energy in the solid phase 
at dimensionless density nr\ = 256, symbols - DMC energy 
(with added tail energy, Eq. (f2)|), solid line - fit E t h + C/VN, 
dashed line - extrapolation of the energy to infinite system 
size E t h- Energy in the inset is measured in units of ti 2 /(Mr 2 ,). 



min{L x , L y }. Here V(r) is the interaction potential and 
92 (f) is the pair distribution function. An approximate 
value of the integral ©: E tail /N = C dd n 3 / 2 /(2\/F) is 
obtained by substituting gi (r) by its average value in the 
bulk g 2 (r) = n. This greatly suppresses finite-size depen- 
dence (for the example shown in the inset of Fig. [1] by a 
factor of 25), while the residual dependence is eliminated 
by a fitting procedure. We note that the dependence 
on the number of particles in this case scales as l/ViV, 
contrary to the law 1/N, usual for fast decaying poten- 
tials. We do the extrapolation of the energy E(N) — 
Edmc + Etaii to the thermodynamic limit value E t h us- 
ing as fitting formula: E(N) = Eth + C/N 1 / 2 , where C is 
a fitting constant. In the inset of Fig. 1 we show an exam- 
ple of the finite size study for the energy. In it we consider 
a density nr^ = 256 and solid phase, where we expect to 
find largest finite-size dependence, due to the oscillatory 
behavior of the pair distribution function. We find that 
our fitting function describes well the finite-size depen- 
dence and we use it for the extrapolation to the thermo- 
dynamic limit. The same procedure is repeated for densi- 
ties nrl = 32; 48; 64; 96; 128; 196; 256; 384; 512; 768; 1024 
and for the gas and solid phases. 

In Fig. Q] we show the results of the ground-state energy 
for gas and solid phases. We find that at small density the 
gas phase is energetically favorable and the solid phase 
is metastable. At larger densities the system crystallizes 
and a triangular lattice is formed. We fit our data points 
with a dependence E/N = ai(nrl) 3 / 2 + a2(nrp) 5 / 4 + 
a 3( nr o) 1 ^ 2 ' where powers 3/2 and 5/4 describe proper 
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FIG. 2: (color online) Static structure factor: left axis - gas 
phase at densities nr\ = 16, 32, 64, 128, 256 (the higher is the 
peak, the larger is the density); right axis - solid phase at 
density nr\ — 256 



asymptotic behavior of potential and kinetic energies 
respectively and power 1/2 is added for a better ac- 
curacy of the fit. The best parameters of the fit are 
oi = 4.536(8); a 2 = 4.38(4); a 3 = 1.2(3) for the gas phase 
and <zi = 4.43(1); o 2 = 4.80(3); a 3 = 2.5(2) for the solid. 
The transition is estimated to happen at the critical den- 
sity nrl = 290(30). We note that in a one-dimensional 
system the role of interactions is enhanced and the tran- 



sition happens at much smaller densities nr w 0.4 14 1. 
Maxwell double tangent construction shows that the re- 
gion of phase coexistence is very small and freezing and 
melting points are indistinguishable within the error bars 
of our calculation. At large densities quantum fluctua- 
tions get suppressed and the energy is dominated by the 
potential energy of particle-particle interactions. The en- 
ergy eventually approaches the limit of a classical crystal. 
The triangular lattice in this limit has potential energy 



/N = 4.446(w ) 3/2 . 



triang 

Knowledge of the equation of state E(n)/N of a ho- 
mogeneous system permits the calculation of the chem- 
ical potential fi(n) = dE/dN and the speed of sound c: 
Ale 2 = nd[i/dn. This information is extremely useful for 
the description of trapped systems with large number of 
particles. As is well known, the local density approxi- 
mation can be used for predictions of the density profile, 
release energy, frequencies of collective modes, etc. which 
can be accessed in experiments [lil ]. 

We estimate the thermodynamic Lindemann ratio 7 = 
yj ((r — r c ) 2 )/aL, (a/, is the lattice length), at the tran- 
sition density as 7 = 0.230(6). Comparing to other two- 
dimensional systems we find that this value is smaller 
than the one of a hard-disks 7 = 0.279, while is similar 
to more long-ranged potentials, like 7 = 0.235(15) for 
Yukawa bosonsand 7 = 0.24(1) for Coulomb bosons [lij. 



In a three-dimensional system the value of 7 at transition 
is generally larger; for example 7 = 0.28 for the Yukawa 
potential 17 1. 

The breaking of translational symmetry takes place in 
the crystal. Periodic modulation appears in the density 
profile n(r) and the pair correlation function gi(r — r') = 
(n(r)n(r')). We have studied the static structure fac- 
tor Sk which is related to <?2 (r) by Fourier transforma- 
tion Sk = 1 + 27r J^°(g2(r) — n)Jo(kr)r dr. We use a 
technique of pure estimators which essentially removes 
any bias from a particular choice of the guiding wave 
function [l8j]. The static structure factor increases from 
zero at k = to unity for large momentum. In the gas 
phase Sk increases smoothly and has a peak around in- 
verse interparticle distance n 1 / 2 (see Fig. As the 
dimensionless density nr 2 , is increased the correlations 
get stronger and the height of the peak is increased. In 
the solid phase static structure factor is discontinuous 
and has a 5-function peak at inverse lattice site length 
al 1 = (V371/2) 1 / 2 = 0.93.. .n 1 ' 2 . 

We can get some insight on the excitation spectrum 
Ek by relating it to the static structure factor through 
the Feynman relation 
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Ek = 



K 2 k 
2MS k ' 



(3) 



This relation is expected to be exact for small momenta 
(i.e. in the phononic regime), while at larger momenta it 
gives an upper bound to the excitation spectrum [l9| . 
Predictions for Ek obtained from DMC data are pre- 
sented in Fig. [3] We find that deviations from low- 
momenta phononic linear behavior appear very soon. 
The excitation spectrum is monotonous for the smallest 
considered density. As the density increases, the roton 
minimum at finite value of the momentum is observed 
(see, also Refs. 0). The roton gap A decreases as the 
density is increased. We expect that our approach gives 
the correct position of the roton minimum, while the real 
value of the roton gap A is overestimated, as it happens 
in DMC calculations for gas helium [20j ]. 

While Bose-Einstein condensation (BEC) is expected 
to happen in a dilute systems, strong interactions be- 
tween particles destroys coherence. In order to study the 
process of decoherence quantitatively we measure con- 
densate density no. In a homogeneous system it can be 
found from the long-range asymptotic \r— r'\ — > 00 of one- 
body density matrix (OBDM) gi(r-r') = (* t (r)*(/)), 
where ^>(r) is the field operator. As g\ (r) is a non-local 
quantity, the DMC output gives for it a mixed estimator, 
which is biased by the choice of the guiding function. 
This bias can be significantly diminished by the extrap- 
olation procedure [18j. We have measured asymptotics of 
OBDM for different sizes of the system and made extrap- 
olation to the thermodynamic limit. In Fig. [4] we plot the 
condensate fraction as a function of the dimensionless 
density in the gas phase. We find that the condensate 
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creases rapidly with the density and is estimated to be 
very small (2%) at the transition point. Our predictions 
can be checked in future experiments with spatially sep- 
arated indirect excitons and cold dipolar Bose gases in 
reduced geometry. 

After completion of this work, we became aware of 
Ref. 2l| in which systems of N=36,90 dipoles are studied 
by finite-temperature Path Integral MC method with a 
result for n c compatible with our prediction. 

Some of the authors (G.E.A. and J.B.) acknowledge 
support by (Spain) Grant No. FIS2005-04181 and Gen- 
eralitat de Catalunya Grant No. 2005SGR-00779. The 
work was partially supported by RFBR. 
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FIG. 3: Upper bound for the excitation spectrum in the gas 
phase. Solid lines - Feynman formula, Eq. <(3j , applied to the 
DMC data for the static structure factor, dashed line - free 
particle limit. Energy is measured in units of h 2 /(Mr ). 
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FIG. 4: Condensate fraction as a function of the dimensionless 
density nr%. 



depletion is large in the range of considered densities. 
This prohibits use of perturbative Bogoliubov approach 
for prediction of the condensate depletion and justifies 
the necessity of a numerical approach to the problem. 
The condensate fraction is rather small at the transition 
density. 

In conclusions, we have investigated the zero- 
temperature quantum phase diagram of a two- 
dimensional dipole system and estimated the critical den- 
sity of the gas to solid phase transition. We studied 
structural properties by observation of the pair distribu- 
tion function, static structure factor, Lindemann ratio. 
Coherence properties are investigated by looking at the 
one-body density matrix. The condensate fraction de- 
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